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Abstract 



We propose a mean-field model for describing the averaged properties 

of a class of stochastic diffusion-limited growth systems. We then show that 

this model exhibits a morphology transition from a dense-branching structure 

with a convex envelope to a dendritic one with an overall concave morphology. 

We have also constructed an order parameter which describes the transition 

quantitatively. The transition is shown to be continuous, which can be verified 

by noting the non-existence of any hysteresis. 
68.70.+w,47.15.Hg,47.20.Hw 
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I. INTRODUCTION 



Diffusion limited growth processes occur in a broad range of interesting systems ranging 
from physics to chemistry and to biology • The common features that we observe in these 
systems are due to the patterns being driven by a generic instability which arises when 
the process of one phase replacing another phase is controlled by a diffusive field [H. This 
instability leads to the breakdown of any simple shape and thereby causes the formation 
of complex interfacial structures. The structures that do form vary from smooth fingers 
to dendritic arbors to disordered, possibly fractal patterns in a manner dependent on the 
details of different individual systems, 

Most of the work in this field has utilized one of two approaches. The first method 
relies on a continuum description of the interface evolution process, leading to a free surface 
problem for partial differential equations governing the diffusive transport. The models that 
one arrives at are purely deterministic; stochastic behavior can of course arise dynamically 
due to the inherent nonlinearities in the interfacial dynamics. A major success of this line 
of reasoning has been the discovery of the "microscopic-solvability" criterion, which enables 
us to understand the selection of a unique finger pattern out of an apparently continuous 
family of available solutions More recently ||, this approach has been used to study 
the global morphology of diffusion-limited solidification. From our perspective, studies of 
solidification using the phase-field method [§J fall into this same category; although quite 
different in computational detail from the sharp interface equations, these methods also 
embody deterministic microscopic models for interfacial evolution. 

The second general scheme available to investigate diffusion-limited growth is the ki- 
netic approach, where the diffusion field is represented by particles executing random walks. 
The boundary condition for the diffusive field and the local dynamics of the interface are 
combined into microscopic rules for these random walkers to stick to the aggregate cluster 
representing the growing "solid" phase. The simplest member of this class of models is dif- 
fusion limited aggregation (DL A) |J , where the sticking rule is simply that the walker will 
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become part of the aggregate upon contact with the aggregate. This simple model is known 
to produce fractal structure and has been extensively investigated 0. However, it is fair 
to say that a full theory of the relationship between the diffusive nature of the controlling 
field (walkers) and the fractal structure of the cluster which DLA generates still eludes us. 
Note that models of this second type are explicitly stochastic; whether this represents an 
important consideration or not is still an open question. 

In order to study more realistic solidification processes, several groups have introduced 
more complicated local kinetics(sticking rules) and in addition have used a finite density of 
random walkers to mimic finite undercooling effects.. Saito and Ueta performed simula- 
tions of a many-walkers model with Ising-like sticking probabilities. Liu and Goldenfeld || 
used a relaxation method in updating the walker distribution function (instead of using walk- 
ers per se) but chose a sticking rule which does not have any well defined thermodynamic 
limit. In a series of recent papers, Shochet et al PJlO| have introduced the "diffusion- 
transition scheme", which they have combined these above two ideas, i.e. directly solving 
for the walker distribution and using an Ising-like sticking probability at the interface. The 
main purpose of this study was to investigate the existence of different patterns as one 
changes the various control parameters, such as the density of the diffusive field at the 
boundary of the system( the undercooling), the chemical potential difference between the 
two phases, (in a spin system analogy, this corresponds to a magnetic field favoring the solid 
phase), and the surface energy (bond energy, in the spin analogy). 

In all the studies mentioned above, a morphology transition from a dense branching 
morphology (DBM) to a dendritic structure has been observed. Locally, the DBM phase 
resembles the ramified structure of a DLA fractal, but at larger length scale, it is densely 
packed, i.e. it has a finite density and the pattern has a well-defined smooth envelop. For the 
dendritic phase, the directions of the dendrites are determined by the anisotropy of either 
the surface energy or the interfacial kinetic coefficient. Due to the probabilistic nature of 
these models, each realization of the pattern looks different; in order to characterize the 
morphology transition, an ensemble average is used to study the statistical properties of 
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the pattern distribution function. In terms of this ensemble averaged pattern, the DBM to 
dendrite transition is connected to the change of the envelop shape from convex to concave. 

Even though such a transition is quite well established in numerical simulations, there 
has been very few theoretical attempts to understand the transition fllTf . Such a study is 
quite crucial if we are to answer questions regarding the role of explicit stochasticity, the 
possibility of morphology selection via a maximum velocity selection principle |12| and in 
fact, the whole concept of a sharp transition between distinct morphologies. 

Our goal here is to propose a new mean field theory (MFT) which describes the ensemble 
averaged behavior of these stochastic models, and more specifically the morphology transi- 
tion. A mean field is again a deterministic set of evolution equations, but these are meant 
to describe the average of the overall probability distribution and not any given realization; 
they should not be confused with the deterministic microscopic models mentioned above. In 
the next section, we describe our previous work on mean field theories for the DLA sticking 
rule and how one can phenomenologically introduce mean field reaction rates which lead to 
a better (i.e. more physically motivated) set of equations. Afterwards, we explicitly discuss 
the issue of the nature of the morphology transition and conclude that in the mean field 
treatment, it is continuous. The final section contains a concluding discussion. 



II. MEAN FIELD EQUATIONS 

In our previous work [Tj|], we have studied the existence and nature of a morphology 
transition for the following set of equations: 

p = uip 1 + a 2 V 2 p) (1) 
u = DV 2 u - p (2) 

where D is the diffusion constant, p is the density of the cluster, u is the density of the 
walkers, with the boundary conditions u(po) = A, p(oo) = 0. Eq. (2) arises due to the con- 
servation of particles and the diffusive nature of the random walkers. Eq.(l) is determined 
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by the local kinetics of the growth process. The above eqs. (1,2) were originally proposed 
to describe ensemble averaged behavior of DLA without the u term in the second equation, 
and with different boundary condition for u(oo) ||14| . The crucial difference between these 
equations and the original Witten-Sander DLA mean field theory is that the phenomenolog- 
ical parameter 7 is taken to be strictly greater than unity; later |T5 ], this type of cutoff in 
the growth rate at small density was derived by proper inclusion of the average effects of the 
multiplicative noise that must be added to the naive reaction kinetics term so as to make 
a complete DLA theory. In that study ||13|| , we did indeed discover that as one lowers the 
undercooling A, there is a morphology transition characterized by the change of the envelop 
shape from convex(DBM) to concave (dendrite). We did not find any discontinuity in the 
velocity slop versus A; this was subsequently verified to be true of the actual transition in 



many- walkers DLA fl6|| . 

It is obvious from eqs. (1,2) that they cannot describe the kinetics of the more complex 
simulations; since the model was supposed to mimic the sticking rules of pure DLA, there is 
no parameter in the mean field model playing the role of the chemical potential difference 
between the two phases. Finding an augmented set of equations is important because the 
aforementioned results on the continuous nature of the transition appear to disagree with 
the findings of the Shochet et al diffusion-transition simulation. Also, a more physical model 
would obviate the need for the parameter 7, since as we shall see, a cutoff in growth at low 
density is an immediate consequence of the actual kinetics of these more realistic simulations. 

We now turn to construction of the new mean field model. One way of understanding 
the crucial feature which we need to include is by recognizing that the local kinetics satisfy 
detailed balance. The local energy functional which governs this balance is determined by 
both the surface energy and the chemical potential difference between the two phases. Thus, 
both the solidification of the liquid and the melting of the solid occur simultaneously with 
their rates dependent on this local energy functional. This leads immediately to the form of 
the simplest MFT as: 
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p = (1 - p)R s (AE, p, s , u) - pR m (AE, p: s , u) 



(3) 



u = DV 2 u - p 



(4) 



The first term on the right hand side of eq. (3) represents the product of the probability 
of the site being empty (1 — p) and the rate of solidification R S (AE, /i s , u); similarly, the 
second represents the product of the probability of the site being occupied by solid p and the 
rate of melting R m (AE, p, s , u). Here AE is the local surface energy and p, s is the chemical 
potential difference between the two phases. 

In the actual kinetic simulation, the surface energy can depend on the detailed geometry 
of the solid near the vicinity of a possible site for melting or solidifying. In the mean field 
description, AE at site i will be taken to depend on pi, which is the average of all nearest- 
neighbor densities at site i. In the continuum limit, p = p + (a 2 /2<f)V 2 p, where a is the 
lattice spacing and d is the dimension. There is no analytic derivation for the explicit form 
of AE(p); however, we know that AE — > oo as p — >• in order to suppress the nucleation 
process inside the liquid phase and AE — >• — oo as p — > 1 to suppress melting inside the solid. 
Also AE = —2B, 0, 2B when p — \, |, | , which reflects the fact that when a site has 1, 2 or 
3 neighbors being occupied, the bond energy gain energy for this site to be occupied is 2B, 
or -2B. We have chosen an expression for AE which satisfies the above requirements: 



We stress that the explicit form of AE(pi) is chosen for convenience, once we ensure that it 
has sensible limits. It will become clear later in the paper that the detailed form of AE is 
unimportant for the qualitative behavior of the system, which in any event is all MFT can 
offer. Once we have an expression for AE, the simplest choices for the transition rates are: 



AE(Pi) = 2B/tg(npi) 



(5) 



(A) R S (AE, pi s ,u) 



u 



R m (AE,p, s ,u) 



1 



(6) 



1 + exp(A£ - fi s ) 



1 + exp(-A£ + fi s ) 



(B) R s (AE,pi s ,u) 



u 



R m (AE,fx s ,u) 



1 



(7) 



u + exp(AE — p, s ) 



1 + u exp(— AE + jig) 
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in accordance with the algorithms used by references and respectively. We can easily 
see that the ratio of the two transition rate are the same for the two algorithms; this is what 
matters for the morphology transition, as we will show later. For the rest of the paper, we 
will use scheme (A). 



III. RESULTS 

We have studied the above equations numerically for both one and two dimensions. 
We first present our results in ID. There are three variables that are important for the 
morphology transition, i. e., B, p s and A. We choose to fix B = 0.7, A = 0.7 and vary the 
chemical potential p s . We find that there is a critical value p* s (B, A), such that when p s > p* 
the dynamics approaches a steady state. That is, starting form any initial condition, the 
system settles into a state in which the solid phase replaces the liquid phase with a constant 
velocity and the profiles of the p and u fields are time independent in the co-moving frame. 
When p s < p*, there is no steady state solution. Instead, the p field will grow up to the 
maximum density p — 1. and the profile of the liquid density is now time-dependent, and 
the width of the u field increases with time. These two phases of eqs. (3,4) are shown in 
figure 1(a), (b). 

Before describing our results for the more interesting 2D case, let us try to understand 
the above results. We look for steady state solution of eq.(3,4), where the interface (or the 
front) moves with certain velocity v, and shape of the field profile does not change with time. 
We transform to the co-moving frame by the change of variable: z = x — vt; the equations 
(3,4) become: 

- Vdp/dz = U{l ~P ] , r (8) 

H/ l + exp(AE-p s ) l + exp(-A£ + /i s ) v ; 

-vdp/dz = Dd 2 u/dz 2 + vdv/dz (9) 
The second equation above can be integrated to give: 

p = A — u — ^-du/dz (10) 



where the boundary conditions p — > 0, u — > A (as z — > oo)has been used. Next, we consider 
the profiles of the p and u fields away from the front, where there is no spatial dependence, 
so: 

u(l - p) - pexpjAE - fi s ) _ 

1 + exp(AE -fi 8 ) ' 1 ' 

p + u = A (12) 

The solution of eq.(ll) is: p = (recall, that AE — > oo as p —> 0) or p = 1 (AE — > — oo) 

or w = i3^ ex P(Ai? — /i s ) (0 < p < 1). The latter nontrivial relation is plotted in figure 2 

together with the straight line p + u = A in the u — p plane. For fi s large enough, there 

are three fixed points A, B and C determined by the eqs. (11,12). Their relative position 

is illustrated in fig. 2, it can be easily seen that fixed points A and C are stable, whereas 

fixed point B is unstable and it separates the attraction basin of the two stable fixed points. 

Fixed point A represents the liquid state where p = 0, u = A, and fixed point C with p ~ A 

and u « A describes the solid state. 

It is well known in non-equilibrium systems fnj, when a stable state invades another 

stable state, the velocity of the front is uniquely determined. One way to see this is to 

substitute the solution of eq. (10) into eq. (8), we have an equation for just u: 

2 m2 r D u + exp(AE - p s ) w(l - A + u) - (A - u) exp(A£ - fi s ) 

Do u/oz + [v — — Aou/oz - 



v 1 + exp(AE — p s ) 1 + exp(A£' — p s ) 

(13) 

The front velocity v is now a parameter in the above equation. By taking into account 
the asymptotic behavior of the u field at z —> ±oo, one can show the above equation is an 
eigenvalue equation for v, i. e., it is only solvable for an unique value of v. 

At some critical plane in (p s , A, B) space, the fixed points B and C disappear together, 
which means there is no steady state growth possible in this regime. In fact, depending on 
the detailed functional form of AE, the steady state solution disappears before the merging 
of B and C when there is no solution for the eigenvalue problem for any v > 0. This is the 
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case for our choice of AE. In the regime where there is no steady state, the solid phase 
flows into the maximum density state with p = 1 (p = 1 is a still a fixed point for eq.(3)), 
and since the liquid phase only supplies a density of it(oo) = A < 1, there is a deficit in 
the u field. The u field tries to compensate for this deficit by getting particles from regions 
increasingly deeper into the liquid side with increasing time and thereby develops a time 
dependent profile. We therefore call this phase of the dynamics the "starved phase", and 
the previous one the "saturated phase". 

Having identified the two phases in ID, we can proceed to study the more important 
2D case. We use the discretized version of eqs. (3,4) using the grid size Ax = Ay = a, 
because any structure which is smaller than the lattice spacing a is unphysical. Let a = 1 
and by writing: p(i,j) = \(p(i + + p(i - I J) + p(j + + p(j - the surface 

energy anisotropy is automatically included in our model. The time step is chosen to be 
small At = 0.01. We find that both of the two phases found in ID have corresponding 
states in 2D. When the chemical potential is large, there is a 2D steady state solution in 
which, in analogy with the "saturated phase" in ID, the solid phase has uniform density 
and the contour line separating the solid and liquid phase has a convex shape (deformed 
from a circle due to anisotropy). This phase can be identified as the DBM phase. The 
2D "starved phase" becomes much more interesting due to the two dimensionality and the 
presence of anisotropy. Because of the surface energy anisotropy, the solid phase has four 
preferred growth directions (45° with respect to the lattice) to grow. As in the ID case, the 
solid phase attempts to grows with maximum density p = 1(> A). However, due to the 
two dimensionality, the deficiency of supply can be compensated by the screening effect, by 
which the growth of solid in the preferred directions screens the growth in other direction, 
and therefore a steady state can be reached. This possibility is directly connected to the 
existence of the Mullins-Serkerka instability. In fact, a stability analysis [TJ[] for our previous 
model shows that there is indeed such instability for the flat interface in the dendritic phase. 
A succession of snaphots of the morphology as we go through the transition by varying the 
chemical potential is shown in fig. 3 . 
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So far, we have shown that eqs. (3,4) have a morphology transition in 2D between two 
phases whose origin can be related to the one dimensional behavior of eqs. (3,4). In order to 
better characterize the transition, we need to have an order parameter which describes the 
different macroscopic nature of the two different morphologies in analogy to what is normally 
done for an equilibrium phase transition. In ||, Shochet et al proposed the front velocity 
as an order parameter; this was based on previous conjectures concerning the selection of 
coexisting patterns Hl2"f . In their simulation, they showed that the dependence of the tip 
velocity on the chemical potential has a discontinuity of slope at the transition point (on 
a log-log plot). We have therefore measured the tip velocity versus the chemical potential 
for our model and the results are shown in figure 4(a), it is quite evident from the plot 
that there is no drastic change at the transition point, in agreement with the results of our 



previous model [jn|. We as yet have no explanation for this discrepancy. 

In general, a useful way to construct an order parameter in non-equilibrium systems is to 
first study the linear stability of the system near the transition. Then the projection of the 
field onto the most unstable mode or the amplitude of the most unstable mode can be used 
as an order parameter. This amplitude equation approach is useful in many non-equilibrium 
systems, including Raleigh Benard convection and Taylor-Couette flow [I7|. However, this 
weakly non-linear methodology does not apply to our case, because dendritic growth is 
highly nonlinear. The initially most unstable mode strongly interacts with other unstable 
modes in the system and the scale of the final pattern is determined by the system size, not 
by the length scale set by the most linearly unstable mode. We are therefore forced to find 
some alternative way to construct an order parameter. 

Given that the very nature of the transition is tied to the change of the spatial pattern, 
it seem natural to us to pick as a measure of the transition a quantity describing how the 
dendritic pattern picks up global correlation as compared to the DBM state. Recall that in 
the dendritic phase, the solid behind the front has a highly non-uniform density; we thus 
define our order parameter as the standard deviation of the solid density: 
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tp = \J< p 2 > - < p > 2 (14) 

with < p 2 >= ^j P fj and < p >= •j P ^ . The integration region uj is the solid dominated 
region defined as where u < A/2. This order parameter ip is essentially the sum of the 
amplitude squares of all the modes. We have plotted the value of ip versus the chemical 
potential in fig. 4(b). It is quite evident from the figure that there is a continuous transition 
at around p,* = 4.8. From the above numerical results, we believe the morphology transition 
(in our model) from the DBM phase to the dendritic phase is continuous. 



To further clarify our results, we would like to comment on the simulations of ref. [10 
which claim to show co-existence of the two distinct morphologies. In that reference, the 
authors performed simulations with the diffusion-transition scheme and confined the growth 
to a channel with the channel direction oriented 45° with respect to the preferred growth 
direction for possible dendrites. In the parameter range where one observes the dendrite 
phase in the open geometry, they obtained a different morphology that is similar to DBM. 
They therefore claimed there exists some range of parameter where both of the two phases are 
"available" , and this coexistence is used to support the idea that transition is discontinuous. 
However, because of the anisotropic nature of the dendritic phase, the boundary condition 
of the channel geometry in |1(| has a very strong effect on the final pattern. It is therefore 
impossible to identify the existence of coexisting phases by using such argument. One way 
to see this is to note that if one continuously changes the angle between the direction of the 
surface energy anisotropy and the boundary, one would see a continuous change of the front 
velocity; this certainly should not be attributed to the existence of a continuous family of 
morphologies. 

To explicitly show how the channel results can be misleading, we have studied our mean 
field model eqs. (3,4) in the channel geometry with the lattice direction chosen to match 
the numerical experiment of ref. [jlO | . We have used channel width W = 50 and channel 
lengths up to L=300. In the parameter regime where DBM exists, we see no change of front 
velocity as shown in fig4(a). This is because DBM phase is a local uncorrelated phase which 
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is not affected by the boundary. When we decrease the chemical potential to the dendrite 
regime, the velocity of the front in the channel diverges from that in the open geometry, 
going below the previous curve. This is due to the strong interaction between the channel 
wall and the dendrites which keep hitting the channel boundary and changing direction. In 
fact, this observation is very useful in identifying the actual morphology transition point; as 
we can see, the transition point determined by this method is the same as determined by 
the order parameter ip. On the other hand, it should not be interpreted as coexistence of 
the two phases, as our model has a continuous transition. 

There is a standard and much more direct way to test whether a transition is continu- 
ous. We can adiabatically change the parameter across the transition point, and determine 
whether there is any hysteresis. We prepare our system in the DBM regime very close to the 
transition point /i s — //*, when we slowly change the chemical potential to lie below //*. The 
pattern immediately starts to generate dendrites at the edge of the previously disordered 
DBM pattern; these grow along the preferred direction, and the forming of the dendrite 
takes over the entire dynamics. A series of plots showing this transition are given in fig. 5. 
The chemical potential /j, s is decreased from /j, s = 5.8 (fig. 5(a))to /j, s = 3.3 (fig. 5(d)). And, 
the reverse process is observed with increasing chemical potential. We thus see no evidence 
for the existence of two attractors with different basins of attraction. 

IV. SUMMARY 

We have proposed a mean field equation to describe the averaged behavior of a class of 
discrete diffusion-limited growth model. Although details of the microscopic kinetics will 
alter the detailed structure of our model, there are universal features which are independent 
of these details. We have shown that as long as there are balancing effects between growth 
and melting, a ID steady state growing solution (the "saturated state") will disappear 
at some critical point; instead, the ID dynamics can be described by the spreading state, 
where the solid grows with the maximum density that is bigger than the supply (the "starved 
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state"). In 2D, these two phases become the DBM phase and the dendrite phase respectively. 
The DBM phase is featureless and the density behind the front is uniform, whereas the 
dendritic phase organizes itself into a correlated structure which then allows higher density 
in parts of the cluster than that could be globally supplied by the liquid state. From these 
properties of the two phases, we have constructed an order parameter, which is the measure 
of uniformity of the solid cluster. The behavior of the order parameter and non-existence of 
hysteresis suggests that the transition is continuous. 

Although we have discounted the purported numerical evidence for the co-existence of 
the DBM and dendritic phases, there still remains the discrepancy between our results on 
the velocity slope and those reported in the simulation studies. Assuming that both our 
results and these kinetic studies withstand further scrutiny, the only remaining possibility 
is that one cannot ignore fluctuation effects even insofar as determining the order of the 
transition. In any case, it is quite clear that mean field approaches can only be used to 
obtain certain types of qualitative information - fluctuation effects are certainly important, 
probably independent of spatial dimensionality, Including these in a complete field-theoretic 
treatment of diffusion-limited growth remains a challenge for future work. 
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FIGURES 



1. The ID profiles of the aggregate and the walker fields with A = 0.7, B = 0.7 for (a) 
the "saturated phase", /i s = 6.8; (b) the "starved phase", /i s = 3.3. 

2. Illustration of the fixed points of the dynamics, fixed points A and C are stable and 
represent the liquid and the solid phases respectively, while fixed point B is unstable 
and separates the attraction basin of the two stable fixed points. 

3. Grey scale plot of the aggregate density field, darker region has lower density. Mor- 
phology changes from DBM (convex) to dendrite (concave) wth fixed A = 0.7, B = 0.7, 
as the chemical potential decreases: (a)// s = 14.3; (b)/i s = 5.8; (c)// s = 3.8; (d)/x s = 2.8. 

4. (a)The front velocity vf versus the chemical potential fi s in a log-log plot. The filled 
circles are data for the open geometry; the empty squares represent the front velocity 
measured in a channel geometry with channel width W = 50. The other parameters 
are fixed as in figure 3: A = 0.7 and B = 0.7. The arrow indicates the morphology 
transition point, (b) The order parameter ip as defined in the text versus the chemical 
potential. 

5. Sequence of grey scale plots of the 2D aggregate density field as the chemical potential 
decreases slowly from above the critical value (DMB regime) /i s — 5.7 (part a) to below 
the critical value (dendrite regime) /i s = 3.3 (part d), part b and c have /i s = 4.9,4.1 
respectively. 
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